1 Carga de paquetes

for (package in c("tidyverse","fpp3", "GGally", "normtest")) {
    if (!require(package, character.only=T, quietly=T)) {
        install.packages(package)
        library(package, character.only=T)
    }
}
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## ── Attaching packages ────────────────────────────────────────────── fpp3 0.3 ──
## ✓ lubridate   1.7.9.2     ✓ feasts      0.1.6  
## ✓ tsibble     0.9.3       ✓ fable       0.2.1  
## ✓ tsibbledata 0.2.0
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## x lubridate::date()   masks base::date()
## x dplyr::filter()     masks stats::filter()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag()        masks stats::lag()
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

2 Funciones auxiliares

# Contraste para los coeficientes
my_t_test <- function (object, ...) 
{
  par <- rbind(t_stat=tidy(object)$statistic, p_value=tidy(object)$p.value)
  colnames(par) <- tidy(object)$term
  if (NCOL(par) > 0) {
    cat("\nt-test:\n")
    coef <- round(par, digits = 4)
    print.default(coef, print.gap = 2)
  }
}

# Gráfico de correlogramas de residuos
my_tsresiduals <- function (data, ...) {
  if (!fabletools::is_mable(data)) {
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  data <- stats::residuals(data)
  if (n_keys(data) > 1) {
    abort("gg_tsresiduals() must be used with a mable containing only one model.")
  }
  gg_tsdisplay(data, !!sym(".resid"), plot_type = "partial",  
    ...)
}

# Validación cruzada anidada

nested_cv <- function(df, h, last_train, string_formula){
  
  nested_errors <- vector()  
  
for (i in seq(last_train, last(df$date), h)){
  train <- df %>% filter(date<=i)
  test <- df %>% filter(date>i)
  
fitted_model <- train %>%
  model(arima = ARIMA(as.formula(string_formula)))

h_forecast = min(dim(test)[1], h)

fc <- fitted_model %>%
  forecast(h=h_forecast)

test_err <- fc %>%
  accuracy(test) %>%
  select(MAPE)
nested_errors <- c(nested_errors, test_err$MAPE)
NewList <- list("errors"=nested_errors, "mean"=mean(nested_errors))
}
return(NewList)}

# Test de autocorrelacines de los residuos
autocorrelation_test_plot <- function(aug, m = 7, dof = 4, h = 5, alpha = 0.05){
vec <- c()
for (i in seq(m, h*m, m)){
 vec <- c(vec,aug %>% features(.resid, ljung_box, lag=i, dof=dof) %>% .$lb_pvalue)
}

autocorr_pvalues_resid <- tibble(
  lag = seq(m, h*m, m), 
  p_value = vec,
  incorelated = p_value >= alpha
)

plot <- autocorr_pvalues_resid %>% 
   ggplot(aes(lag, p_value, color = incorelated)) + 
  geom_point() + 
  geom_hline(aes(yintercept = alpha), linetype="dashed", color = "indianred2")


newList <- list("values" = autocorr_pvalues_resid, "plot" = plot)
return(newList)
}

3 División train y test

Antes de empezar a hacer el estudio hay que distinguir la parte de los datos que se utilizará para construir la fórmula (para modelizar) de aquella que se utilizará para validar los resultados y cuyos datos, no deben ser utilizados.

Convertimos la serie en un objeto tsibble, para trabajar de manera cómoda. Analizamos visualmente la serie y las marcas de tiempo para saber cómo dividir los datos.

demandaGas <- read_csv('DemandaGas.csv')
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   date = col_date(format = ""),
##   value = col_double()
## )
demandaGas = as_tsibble(demandaGas, index = date)

demandaGas %>%
  autoplot(value) +
    labs(title = "Daily Gas Demand") +
    xlab("Year") + ylab("value")

Dado que tenemos datos mensuales y parece existir a simple vista un patrón estacional anual, reservamos los dos últimos años para test y el resto para train.

demandaGas_train <- demandaGas %>% filter_index(. ~ "1999-06-30")
demandaGas_test <- demandaGas %>% filter_index("1999-07-01" ~ .)

4 Fase de identificación

Se comienza la fase de identificación donde se deciden las transformaciones a realizar y el propio ajuste de la serie. Aunque veremos que cuando se realiza la selección de los parámetros del modelo ARIMA, la iteración con la estimación y el contraste son constantes hasta converger a un modelo válido.

4.1 Estacionariedad

Antes de analizar la serie mediante test estadísticos, analicemos gráficamente más en detalle, para confirmar la estacionalidad detectada.

demandaGas_train %>%
gg_season(value, labels = "right")

Parece claro mediante este gráfico que hay una componente estacional anual, dado que la forma de las series anuales son iguales. Además, se observa que el valor aumenta cada año por lo que parece que tampoco es estable en media. También, es altamente probable que exista estacionalidad semanal, dado que la serie es diaria y los consumos de gas variarán en los días laborables y no laborables.

Dadas las características de la serie, también podemos calcular la descomposición STL para corroborar que existe tendencia y estacionalidad.

demandaGas_train %>%
   model(seats = feasts:::STL(value)) %>%
  components() %>%
  autoplot()

Se confirman visualmente las estacionalidades supuestas. Además de la tendencia creciente. Por tanto, parece que la serie no va a ser estacionaria en media. Como era de esperar.

Estacionariedad en varianza

Analicemos qué transformación Box-Cox puede aplicarse para estabilizar la varianza.

lambda <- demandaGas_train %>%
features(value, features = guerrero) %>% pull(lambda_guerrero)
lambda
## [1] 0.6617215

Dado que es un valor próximo a 1, no realizaremos ningún cambio a priori. Si en el proceso de estimaicón no conseguimos un modelo adecuado, volveríamos a este punto para valorar de nuevo la transformación.

Estacionariedad en media

Es evidente que se necesita realizar al menos una diferencia para estabilizar la media. Esto se puede comprobar mediante una prueba de raices unitarias.

demandaGas_train %>%
  features(value, unitroot_kpss)

El p-valor es menor que 0.05, lo que indica que la hipótesis nula es rechazada. Es decir, los datos no son estacionarios. Se confirma formalmente que los datos no son estacionarios.

demandaGas_train %>%
  features(difference(value, 1), unitroot_kpss)

La serie ya pasa el test, por tanto se puede concluir que hay que realizar una diferencia regular.

demandaGas_train %>% autoplot(difference(value, 1))

Debido a la estacionalidad de los datos, quizá deba realizarse una diferencia estacional. Para evaluarlo podemos recurrir a la función unitroot_nsdiffs(), esta evalúa el estadístico fuerza estacional \(F_S\) y sugiere una diferencia estacional si esta es mayor que 0.64.

demandaGas_train %>%
  mutate(turnover = difference(value,1)) %>%
  features(turnover, unitroot_nsdiffs)

Efectivamente, el método sugiere una diferencia estacional como habíamos sospechado.

4.2 Determinación del modelo

En los pasos previos hemos detectado la necesidad de realizar una diferencia regular y otra estacional. Por tanto, comenzaremos con el ajuste de un modelo SARIMA (0,1,0)x(0,1,0)\(_{7}\), sobre la serie.

fit <- demandaGas_train %>%
  model(arima = ARIMA(value ~ pdq(0,1,0) + PDQ(0,1,0, period = 7)))
fit %>% my_tsresiduals(lag_max = 36)

El nuevo gráfico acf corresponde a un AR(7). Se aprecia estructura multiplicativa. Es decir, un AR(1) estacional de periodo 7.

5 Fase de estimación y contraste

Comencemos ajustando la parte estacional mediante un SARIMA (0,1,0)x(1,1,0)\(_{7}\)

fit2 <- demandaGas_train %>%
  model(arima = ARIMA(value ~ pdq(0,1,0) + PDQ(1,1,0, period=7)))
fit2 %>% my_tsresiduals(lag_max =36)

report(fit2)
## Series: value 
## Model: ARIMA(0,1,0)(1,1,0)[7] 
## 
## Coefficients:
##          sar1
##       -0.3450
## s.e.   0.0284
## 
## sigma^2 estimated as 90082:  log likelihood=-7742.84
## AIC=15489.67   AICc=15489.68   BIC=15499.65
my_t_test(fit2)
## 
## t-test:
##              sar1
## t_stat   -12.1542
## p_value    0.0000

Hemos comprobado que el parámetro de este modelo es significativo y no existen problemas de invertivilidad.

Siguen existiendo autocorrelaciones fuertes de orden 7. Por tanto, ajustaremos un MA(7), al resto del modelo. Por tanto, SARIMA (0,1,0)x(1,1,1)\(_{7}\)

fit3 <- demandaGas_train %>%
  model(arima = ARIMA(value ~ pdq(0,1,0) + PDQ(1,1,1, period=7)))
fit3 %>% my_tsresiduals(lag_max =36)

Parece observarse una estructura AR(\(\inf\)) en la PACF. Por tanto, observamos el orden de ese MA finito en la ACF, parece ser 2. Ajustaremos en el siguiente paso un MA(2) regular.

Comprobemos la calidad estadística de este nuevo modelo.

report(fit3)
## Series: value 
## Model: ARIMA(0,1,0)(1,1,1)[7] 
## 
## Coefficients:
##         sar1     sma1
##       0.1938  -0.9670
## s.e.  0.0310   0.0085
## 
## sigma^2 estimated as 62986:  log likelihood=-7555.68
## AIC=15117.35   AICc=15117.37   BIC=15132.33
my_t_test(fit3)
## 
## t-test:
##            sar1       sma1
## t_stat   6.2592  -114.2869
## p_value  0.0000     0.0000

Vemos que de nuevo todos los parámetros son significativos. Otra forma de evaluar si hay problemas con la invertibilidad o estacionariedad es comprobar si el valor estimado del parámetro \(\pm\) el error estándar, en valor absoluto, continene al 1. Por tanto, en este caso no hay problemas con la significatividad de los parámetros ni con la condición de invertibilidad.

Añadimos el MA(2) regular al resto del modelo. Por tanto, SARIMA (0,1,2)x(1,1,1)\(_{7}\)

fit4 <- demandaGas_train %>%
  model(arima = ARIMA(value ~ pdq(0,1,2) + PDQ(1,1,1)))
fit4 %>% my_tsresiduals(lag_max =36)

report(fit4)
## Series: value 
## Model: ARIMA(0,1,2)(1,1,1)[7] 
## 
## Coefficients:
##          ma1      ma2    sar1     sma1
##       0.0905  -0.1875  0.1904  -0.9683
## s.e.  0.0299   0.0306  0.0310   0.0084
## 
## sigma^2 estimated as 60596:  log likelihood=-7533.84
## AIC=15077.68   AICc=15077.73   BIC=15102.64
my_t_test(fit4)
## 
## t-test:
##             ma1      ma2    sar1       sma1
## t_stat   3.0245  -6.1282  6.1446  -114.7039
## p_value  0.0025   0.0000  0.0000     0.0000

No es posible identificar más estructura a partir del ACF y PACF, no hay problemas de estacionariedad o invertibilidad. Además, todos los coeficientes son significativos. Aunque algunas de las autocorrelaciones se siguen saliendo de las bandas.

NOTA: Hay un aspecto que pasa inadvertido por la propia construcción del paquete empleado. No se ha evaluado la posibilidad de incluir una constante en el modelo. Esto es, porque al no pasar el t-test en ningún modelo evaluado directamente no se calcula. Para incluir esta constate y poder comprobarlo, simplemente incluímos “+1” en la fórmula del ARIMA.

fit5 <- demandaGas_train %>%
  model(arima = ARIMA(value ~ pdq(0,1,2) + PDQ(1,1,1) + 1))
## Warning: Model specification induces a quadratic or higher order polynomial trend. 
## This is generally discouraged, consider removing the constant or reducing the number of differences.
fit5 %>% my_tsresiduals(lag_max =36)

report(fit5)
## Series: value 
## Model: ARIMA(0,1,2)(1,1,1)[7] w/ poly 
## 
## Coefficients:
##          ma1      ma2    sar1     sma1  constant
##       0.0905  -0.1875  0.1904  -0.9683   -0.0123
## s.e.  0.0299   0.0306  0.0310   0.0084    0.2724
## 
## sigma^2 estimated as 60652:  log likelihood=-7533.84
## AIC=15079.68   AICc=15079.75   BIC=15109.62
my_t_test(fit5)
## 
## t-test:
##             ma1      ma2    sar1       sma1  constant
## t_stat   3.0247  -6.1281  6.1447  -114.6897   -0.0451
## p_value  0.0025   0.0000  0.0000     0.0000    0.9641

Se aprecia que no podemos rechazar la hipótesis de que la constante sea nula, por tanto nos quedamos con el modelo anterior sin constante. Procedemos al análisis de los residuos del modelo SARIMA (0,1,2)x(1,1,1)\(_{7}\).

5.1 Diagnosis de residuos

Para completar la fase de contraste evluamos los residuos. Veamos primero el histograma.

aug <-fit4 %>% augment()

# Histogram
aug %>%
  ggplot(aes(x = .resid)) +
  geom_histogram(bins = 50) +
  ggtitle("Histogram of residuals")

Test media = 0

Parece que la media es cero, pero distan bastante de la normalidad a priori. Procedemos a realizar un contraste t-student para contrastar si la media es 0.

# Student's t-Test for mean=0
t.test(aug$.resid)
## 
##  One Sample t-test
## 
## data:  aug$.resid
## t = -0.020738, df = 1094, p-value = 0.9835
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -14.67620  14.36923
## sample estimates:
##  mean of x 
## -0.1534891

Como p-value > 0.05 no podemos rechazar la hipótesis de que la muestra tiene media 0.

Test autocorrelaciones

A continuación comprobamos que los residuos están incorrelados.

# Ljung-Box autocorrelation
autocorrelation_test_plot(aug, m = 7, dof = 4, h = 5, alpha = 0.05)
## $values
## # A tibble: 5 x 3
##     lag   p_value incorelated
##   <dbl>     <dbl> <lgl>      
## 1     7 0.0812    TRUE       
## 2    14 0.000231  FALSE      
## 3    21 0.0000496 FALSE      
## 4    28 0.000654  FALSE      
## 5    35 0.000215  FALSE      
## 
## $plot

Los p-valores son menores de 0.05 en la mayoría de las autocorrelaciones. Por lo tanto, NO podemos concluir que los residuos no están correlados.

Test homocedasticidad

Para contrastar la heterocedasticidad se puede utilizamos una regersión media-dispersión. Calculando los grupos de manera anual, dado que esa es la estacionalidad de la serie.

log_log <- aug %>% as_tibble() %>% 
  group_by(week(date)) %>% 
  summarize(mean_resid = log(mean(.resid+1)), std_resid = log(sd(.resid+1))) 
## `summarise()` ungrouping output (override with `.groups` argument)
  summary(lm(std_resid~mean_resid, log_log))
## 
## Call:
## lm(formula = std_resid ~ mean_resid, data = log_log)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7367 -0.3183 -0.1087  0.2913  1.2887 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.00560    0.19332  25.892   <2e-16 ***
## mean_resid   0.04686    0.05966   0.785     0.44    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5273 on 24 degrees of freedom
##   (27 observations deleted due to missingness)
## Multiple R-squared:  0.02506,    Adjusted R-squared:  -0.01556 
## F-statistic: 0.617 on 1 and 24 DF,  p-value: 0.4399

Como el p-valor de log(media) es grande, entonces no puede considerarse distinto de 0 y por tanto los residuos son homocedásticos.

Test de normalidad

Para contrastar la normalidad realizamos el test Jarque-Bera para evaluar la normalidad.

# Jarque Bera test
jb.norm.test(na.omit(aug$.resid))
## 
##  Jarque-Bera test for normality
## 
## data:  na.omit(aug$.resid)
## JB = 7171, p-value < 2.2e-16

Como el p-valor es menor que el nivel de significancia 0.05, se rechaza la hipótsis nula de normalidad de los resíduos como habíamos supuesto por el histograma.

6 Fase de predicción

Finalmente evaluamos la capacidad predictiva del modelo. Calculamos las predicciones en el conjunto de test, calculamos sus errores y los comparamos con los residuos.

# Residual accuracy
resids <- fit4 %>% 
  accuracy() %>% 
  select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>% 
  mutate(Evaluation='Training') 

# Forecasting
fc <- fit4 %>%
  forecast(h=7) 

test_err <- fc %>% 
  accuracy(demandaGas) %>% 
  select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>% 
  mutate(Evaluation='Test')

# Show errors together
bind_rows(test_err, resids) %>% select(Evaluation, everything())

Además de evaluar los errores con estas medidas globales, podemos evaluar los errores cometidos año a año visualmente. Así, podemos detectar outliers o efectos de calendario no detectados. Pudiendo así corregirlos con modelos más complejos mediante variables exógenas.

aug %>%  ggplot() +
  geom_line(aes(x = date, y = .fitted), color="navy") +
  geom_line(aes(x = date, y = value), color="gray24") +
  # geom_line(data=demandaGas_forecast, aes(x = index, y = value), color="red4") +
  ggtitle("SARIMA train fitted values") +
  xlab('Dates') +
  ylab('Demand') + facet_wrap(vars(year(date)), scales = 'free')

aug %>% filter(year(date) == 1998) %>% 
  ggplot() +
  geom_line(aes(x = date, y = .fitted), color="navy") +
  geom_line(aes(x = date, y = value), color="gray24") +
  # geom_line(data=demandaGas_forecast, aes(x = index, y = value), color="red4") +
  ggtitle("SARIMA train fitted values") +
  xlab('Dates') +
  ylab('Demand') + facet_wrap(vars(month(date)), scales = 'free')

También podemos dibujar las predicciones con los intervalos de confianza (aunque en este caso no sean representativos, dada la no normalidad de los residuos). Evaluamos así también la capacidad de generalización del modelo, comparando las predicciones con el conjunto de test.

h <- dim(demandaGas_test)[1]

demanda_plot <- demandaGas %>% filter(date>last(demandaGas_train$date)-14 & date< last(demandaGas_train$date) + 7 )

fit4 %>%
  forecast(h=7) %>%
  autoplot(demanda_plot)

Como sabemos que los modelos SARIMA no predicen bien más alla del orden más alto disponible, en este caso 7, los errores en el conjunto de test más allá de 7 datos se dipararían. Es un claro ejemplo donde puede realizarse nested cross-validation.

A continuación, se muestra el error medio en ventanas deslizantes de tamaño 7 sin solapamiento.

nested_av_errors <- nested_cv(demandaGas, 7, last(demandaGas_train$date), "value ~ pdq(0,1,2) + PDQ(1,1,1)")
nested_av_errors
## $errors
##  [1]  4.9432786  4.5946162  0.4544413  4.7910943  4.6278381  8.2019157
##  [7]  8.7199033 10.9214904 14.4622008  8.6540625  1.8901096  4.6460032
## [13]  2.1681745  1.5557131  3.7169177  2.1928115  2.7725992  3.4787309
## [19]  4.2783305  1.4240245  1.5669736  2.1543376  6.0760687 10.9050842
## [25]  1.7995028 52.4119725 15.5534924  9.9182339  2.3166850  2.5544861
## [31]  4.6191310  6.3948084  2.1454387  2.1260729  2.7343799  2.1228445
## [37]  2.5592099  1.6424454  4.1090883  1.2036735  2.6514637  2.3181956
## [43] 24.4984560  9.5457710  2.0803222  3.2394402  1.0640503  3.3559921
## [49]  3.1263909  2.9427594  1.7663514  1.3239572  1.1677175
## 
## $mean
## [1] 5.518661

Con este error, podemos concluir que una estimación más realista del error de nuestro modelo es de un 5.5%.